!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine objective_shapeener_dx(co,kon,ipkon,lakon,ne,
     &  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
     &  ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
     &  iprestr,iperturb,iout,vold,nmethod,
     &  veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
     &  xstateini,xstiff,xstate,npmat_,matname,mi,ielas,icmd,
     &  ncmat_,nstate_,stiini,vini,ener,enerini,istep,iinc,
     &  springarea,reltime,calcul_qa,iener,ikin,ne0,thicke,
     &  emeini,pslavsurf,pmastsurf,mortar,clearini,nea,neb,ielprop,
     &  prop,distmin,ndesi,nodedesi,nobject,g0,dgdx,
     &  iobject,sti,xener,istartdesi,ialdesi,xdesi,idesvar)
!
!     calculates the internal energy and its derivative w.r.t. the coordinates
!     of the mesh
!
      implicit none
!
      integer cauchy
!
      character*8 lakon(*),lakonl
      character*80 amat,matname(*)
!
      integer kon(*),konl(26),nea,neb,mi(*),mint2d,nopes,
     &  nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),
     &  ielorien(mi(3),*),ntmat_,ipkon(*),ne0,iflag,null,
     &  istep,iinc,mt,ne,mattyp,ithermal(2),iprestr,i,ii,j,k,m1,m2,jj,
     &  i1,m3,kk,iener,indexe,nope,norien,iperturb(*),iout,
     &  icmd,ihyper,nmethod,kode,imat,mint3d,iorien,ielas,
     &  istiff,ncmat_,nstate_,ikin,ilayer,nlayer,ki,kl,ielprop(*),
     &  nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_qa,
     &  nopered,mortar,jfaces,iloc,igauss,idesvar,ndesi,nodedesi(*),
     &  nobject,iobject,iactive,node,istartdesi(*),ialdesi(*),ij
!
      real*8 co(3,*),shp(4,26),stiini(6,mi(1),*),xener(*),
     &  stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*),
     &  elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),xs2(3,7),
     &  alcon(0:6,ntmat_,*),vini(0:mi(2),*),thickness,
     &  alzero(*),orab(7,*),elas(21),rho,qa(3),
     &  fnl(3,10),beta(6),q(0:mi(2),26),xl2(3,8),
     &  vkl(0:3,3),t0(*),t1(*),prestr(6,mi(1),*),
     &  ckl(3,3),vold(0:mi(2),*),eloc(9),veold(0:mi(2),*),
     &  springarea(2,*),elconloc(21),eth(6),xkl(3,3),voldl(0:mi(2),26),
     &  xikl(3,3),ener(mi(1),*),emec(6),enerini(mi(1),*),
     &  emec0(6),veoldl(0:mi(2),26),xsj2(3),shp2(7,8),
     &  e,un,al,um,am1,xi,et,ze,tt,exx,eyy,ezz,exy,exz,eyz,
     &  xsj,vj,t0l,t1l,dtime,weight,pgauss(3),vij,time,ttime,
     &  plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
     &  xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802),
     &  vokl(3,3),xstateini(nstate_,mi(1),*),vikl(3,3),
     &  gs(8,4),a,reltime,tlayer(4),dlayer(4),xlayer(mi(3),4),
     &  thicke(mi(3),*),emeini(6,mi(1),*),clearini(3,9,*),
     &  pslavsurf(3,*),pmastsurf(6,*),distmin,xenerel,g0(iobject),
     &  dgdx(ndesi,nobject),sti(6,mi(1),*),xdesi(3,*)
!
      include "gauss.f"
!
      iflag=3
      null=0
!
      mt=mi(2)+1
!
!     -------------------------------------------------------------
!     Initialisation of the loop for the variation of 
!     the internal forces
!     -------------------------------------------------------------
!
!     
!     Loop over all elements in thread
         do ij=nea,neb
            if(idesvar.gt.0) then
               i=ialdesi(ij)
            else
               i=ij
            endif
!     
!     initialisation of xenerel and 
!     
            xenerel=0.d0
            lakonl=lakon(i)
!     
            if(ipkon(i).lt.0) cycle
!     
!     no 3D-fluid elements
!     
            if(lakonl(1:1).eq.'F') cycle
            if(lakonl(1:7).eq.'DCOUP3D') cycle
!     
            if(lakonl(7:8).ne.'LC') then
!     
               imat=ielmat(1,i)
               amat=matname(imat)
               if(norien.gt.0) then
                  iorien=ielorien(1,i)
               else
                  iorien=0
               endif
!     
               if(nelcon(1,imat).lt.0) then
                  ihyper=1
               else
                  ihyper=0
               endif
            else
!     
!     composite materials
!     
               mint2d=4
               nopes=8
!     determining the number of layers
!     
               nlayer=0
               do k=1,mi(3)
                  if(ielmat(k,i).ne.0) then
                     nlayer=nlayer+1
                  endif
               enddo
!     
!     determining the layer thickness and global thickness
!     at the shell integration points
!     
               iflag=1
               indexe=ipkon(i)
               do kk=1,mint2d
                  xi=gauss3d2(1,kk)
                  et=gauss3d2(2,kk)
                  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
                  tlayer(kk)=0.d0
                  do k=1,nlayer
                     thickness=0.d0
                     do j=1,nopes
                        thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
                     enddo
                     tlayer(kk)=tlayer(kk)+thickness
                     xlayer(k,kk)=thickness
                  enddo
               enddo
               iflag=3
!     
               ilayer=0
               do k=1,4
                  dlayer(k)=0.d0
               enddo
!     
            endif
!     
            indexe=ipkon(i)
c     Bernhardi start
            if(lakonl(1:5).eq.'C3D8I') then
               nope=11
            elseif(lakonl(4:5).eq.'20') then
c     Bernhardi end
               nope=20
            elseif(lakonl(4:4).eq.'2') then
               nope=26
            elseif(lakonl(4:4).eq.'8') then
               nope=8
            elseif(lakonl(4:5).eq.'10') then
               nope=10
            elseif(lakonl(4:5).eq.'14') then
               nope=14
            elseif(lakonl(4:4).eq.'4') then
               nope=4
            elseif(lakonl(4:5).eq.'15') then
               nope=15
            elseif(lakonl(4:4).eq.'6') then
               nope=6
            elseif((lakonl(1:1).eq.'E').and.(lakonl(7:7).ne.'F')) then
!     
!     spring elements, contact spring elements and
!     dashpot elements
!     
               if(lakonl(7:7).eq.'C') then
!     
!     contact spring elements
!     
                  if(mortar.eq.1) then
!     
!     face-to-face penalty
!     
                     nope=kon(ipkon(i))
                  elseif(mortar.eq.0) then
!     
!     node-to-face penalty
!     
                     nope=ichar(lakonl(8:8))-47
                     konl(nope+1)=kon(indexe+nope+1)
                  endif
               else
!     
!     genuine spring elements and dashpot elements
!     
                  nope=ichar(lakonl(8:8))-47
               endif
            else
               cycle
            endif
!     
!     check whether a design variable belongs to the element    
!     
            if(idesvar.gt.0) then 
               do ii=1,nope
                  node=kon(indexe+ii)
                  if(node.eq.nodedesi(idesvar)) then
                     iactive=ii
                     exit
                  endif
               enddo
            endif
!     
            if(lakonl(4:5).eq.'8R') then
               mint3d=1
            elseif(lakonl(4:7).eq.'20RB') then
               if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
                  mint3d=50
               else
                  call beamintscheme(lakonl,mint3d,ielprop(i),prop,
     &                 null,xi,et,ze,weight)
               endif
            elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'26R').or.
     &              (lakonl(4:6).eq.'20R')) then
               if(lakonl(7:8).eq.'LC') then
                  mint3d=8*nlayer
               else
                  mint3d=8
               endif
            elseif(lakonl(4:4).eq.'2') then
               mint3d=27
            elseif((lakonl(4:5).eq.'10').or.(lakonl(4:5).eq.'14')) then
               mint3d=4
            elseif(lakonl(4:4).eq.'4') then
               mint3d=1
            elseif(lakonl(4:5).eq.'15') then
               mint3d=9
            elseif(lakonl(4:4).eq.'6') then
               mint3d=2
            elseif(lakonl(1:1).eq.'E') then
               mint3d=0
            endif
            
            do j=1,nope
               konl(j)=kon(indexe+j)
               do k=1,3
                  xl(k,j)=co(k,konl(j))
                  vl(k,j)=vold(k,konl(j))
                  voldl(k,j)=vold(k,konl(j))
               enddo
            enddo
!     
!     computation of the objective function and the derivate 
!     
            if(idesvar.gt.0) then 
               do j=1,3
                  xl(j,iactive)=
     &                 xl(j,iactive)+xdesi(j,idesvar)
               enddo
            endif
!     
!     calculating the forces for the contact elements
!     
            if(mint3d.eq.0) then
!     
!     "normal" spring and dashpot elements
!     
               kode=nelcon(1,imat)
               if(lakonl(7:7).eq.'A') then
                  t0l=0.d0
                  t1l=0.d0
                  if(ithermal(1).eq.1) then
                     t0l=(t0(konl(1))+t0(konl(2)))/2.d0
                     t1l=(t1(konl(1))+t1(konl(2)))/2.d0
                  elseif(ithermal(1).ge.2) then
                     t0l=(t0(konl(1))+t0(konl(2)))/2.d0
                     t1l=(vold(0,konl(1))+vold(0,konl(2)))/2.d0
                  endif
               endif
!     
!     spring elements (including contact springs)
!     
               if(lakonl(2:2).eq.'S') then
                  if((lakonl(7:7).eq.'A').or.((mortar.eq.0).and.
     &                 ((nmethod.ne.1).or.(iperturb(1).ge.2).or.
     &                  (iout.ne.-1))))then
                     call springforc_n2f(xl,konl,vl,imat,elcon,nelcon,
     &                    elas,fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,
     &                    elconloc,plicon,nplicon,npmat_,ener(1,i),
     &                    iener,stx(1,1,i),mi,springarea(1,konl(nope+1))
     &                    ,nmethod,ne0,nstate_,xstateini,xstate,reltime,
     &                    ielas,ener(1,i+ne))
                  elseif((mortar.eq.1).and.
     &                    ((nmethod.ne.1).or.(iperturb(1).ge.2).or.
     &                     (iout.ne.-1)))then
                     iloc=kon(indexe+nope+1)
                     jfaces=kon(indexe+nope+2)
                     igauss=kon(indexe+nope+1)
                     call springforc_f2f(xl,vl,imat,elcon,nelcon,elas,
     &                    fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,
     &                    elconloc,plicon,nplicon,npmat_,ener(1,i),
     &                    iener,stx(1,1,i),mi,springarea(1,iloc),
     &                    nmethod,ne0,nstate_,xstateini,xstate,reltime,
     &                    ielas,iloc,jfaces,igauss,pslavsurf,pmastsurf,
     &                    clearini,ener(1,i+ne))
                  endif
!     
               endif
c            elseif(ikin.eq.1) then
c               do j=1,nope
c                  do k=1,3
c                     veoldl(k,j)=veold(k,konl(j))
c                  enddo
c               enddo            
            endif
!     
            do jj=1,mint3d
!     
               if(lakonl(4:5).eq.'8R') then
                  xi=gauss3d1(1,jj)
                  et=gauss3d1(2,jj)
                  ze=gauss3d1(3,jj)
                  weight=weight3d1(jj)
               elseif(lakonl(4:7).eq.'20RB') then
                  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
                     xi=gauss3d13(1,jj)
                     et=gauss3d13(2,jj)
                     ze=gauss3d13(3,jj)
                     weight=weight3d13(jj)
                  else
                     call beamintscheme(lakonl,mint3d,ielprop(i),prop,
     &                    jj,xi,et,ze,weight)
                  endif
               elseif((lakonl(4:4).eq.'8').or.
     &                 (lakonl(4:6).eq.'20R').or.(lakonl(4:6).eq.'26R'))
     &                 then
                  if(lakonl(7:8).ne.'LC') then
                     xi=gauss3d2(1,jj)
                     et=gauss3d2(2,jj)
                     ze=gauss3d2(3,jj)
                     weight=weight3d2(jj)
                  else
                     kl=mod(jj,8)
                     if(kl.eq.0) kl=8
!     
                     xi=gauss3d2(1,kl)
                     et=gauss3d2(2,kl)
                     ze=gauss3d2(3,kl)
                     weight=weight3d2(kl)
!     
                     ki=mod(jj,4)
                     if(ki.eq.0) ki=4
!     
                     if(kl.eq.1) then
                        ilayer=ilayer+1
                        if(ilayer.gt.1) then
                           do k=1,4
                              dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
                           enddo
                        endif
                     endif
                     ze=2.d0*(dlayer(ki)+(ze+1.d0)/
     &                    2.d0*xlayer(ilayer,ki))/tlayer(ki)-1.d0
                     weight=weight*xlayer(ilayer,ki)/tlayer(ki)
!     
!     material and orientation
!     
                     imat=ielmat(ilayer,i)
                     amat=matname(imat)
                     if(norien.gt.0) then
                        iorien=ielorien(ilayer,i)
                     else
                        iorien=0
                     endif
!     
                     if(nelcon(1,imat).lt.0) then
                        ihyper=1
                     else
                        ihyper=0
                     endif
                  endif
               elseif(lakonl(4:4).eq.'2') then
                  xi=gauss3d3(1,jj)
                  et=gauss3d3(2,jj)
                  ze=gauss3d3(3,jj)
                  weight=weight3d3(jj)
               elseif((lakonl(4:5).eq.'10').or.(lakonl(4:5).eq.'14'))
     &            then
                  xi=gauss3d5(1,jj)
                  et=gauss3d5(2,jj)
                  ze=gauss3d5(3,jj)
                  weight=weight3d5(jj)
               elseif(lakonl(4:4).eq.'4') then
                  xi=gauss3d4(1,jj)
                  et=gauss3d4(2,jj)
                  ze=gauss3d4(3,jj)
                  weight=weight3d4(jj)
               elseif(lakonl(4:5).eq.'15') then
                  xi=gauss3d8(1,jj)
                  et=gauss3d8(2,jj)
                  ze=gauss3d8(3,jj)
                  weight=weight3d8(jj)
               elseif(lakonl(4:4).eq.'6') then
                  xi=gauss3d7(1,jj)
                  et=gauss3d7(2,jj)
                  ze=gauss3d7(3,jj)
                  weight=weight3d7(jj)
               endif
!     
c     Bernhardi start
               if(lakonl(1:5).eq.'C3D8R') then
                  call shape8hr(xl,xsj,shp,gs,a)
               elseif(lakonl(1:5).eq.'C3D8I') then
                  call shape8hu(xi,et,ze,xl,xsj,shp,iflag)
               elseif(nope.eq.20) then
c     Bernhardi end
                  if(lakonl(7:7).eq.'A') then
                     call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
                  elseif((lakonl(7:7).eq.'E').or.
     &                    (lakonl(7:7).eq.'S')) then
                     call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
                  else
                     call shape20h(xi,et,ze,xl,xsj,shp,iflag)
                  endif
               elseif(nope.eq.26) then
                  call shape26h(xi,et,ze,xl,xsj,shp,iflag,konl)
               elseif(nope.eq.8) then
                  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
               elseif(nope.eq.10) then
                  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
               elseif(nope.eq.14) then
                  call shape14tet(xi,et,ze,xl,xsj,shp,iflag,konl)
               elseif(nope.eq.4) then
                  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
               elseif(nope.eq.15) then
                  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
               else
                  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
               endif
!     
!     vkl(m2,m3) contains the derivative of the m2-
!     component of the displacement with respect to
!     direction m3
!     
               do m2=1,3
                  do m3=1,3
                     vkl(m2,m3)=0.d0
                  enddo
               enddo
!     
               do m1=1,nope
                  do m2=1,3
                     do m3=1,3
                        vkl(m2,m3)=vkl(m2,m3)+shp(m3,m1)*vl(m2,m1)
                     enddo
                  enddo
               enddo
!     
!     for frequency analysis or buckling with preload the
!     strains are calculated with respect to the deformed
!     configuration
!     for a linear iteration within a nonlinear increment:
!     the tangent matrix is calculated at strain at the end
!     of the previous increment
!     
               if((iperturb(1).eq.1).or.(iperturb(1).eq.-1))then
                  do m2=1,3
                     do m3=1,3
                        vokl(m2,m3)=0.d0
                     enddo
                  enddo
!     
                  do m1=1,nope
                     do m2=1,3
                        do m3=1,3
                           vokl(m2,m3)=vokl(m2,m3)+
     &                          shp(m3,m1)*voldl(m2,m1)
                        enddo
                     enddo
                  enddo
               endif
!     
               kode=nelcon(1,imat)
!     
!     calculating the strain
!     
!     attention! exy,exz and eyz are engineering strains!
!     
               exx=vkl(1,1)
               eyy=vkl(2,2)
               ezz=vkl(3,3)
               exy=vkl(1,2)+vkl(2,1)
               exz=vkl(1,3)+vkl(3,1)
               eyz=vkl(2,3)+vkl(3,2)
!     
               if(iperturb(2).eq.1) then
!     
!     Lagrangian strain
!     
                  exx=exx+(vkl(1,1)**2+vkl(2,1)**2+vkl(3,1)**2)/2.d0
                  eyy=eyy+(vkl(1,2)**2+vkl(2,2)**2+vkl(3,2)**2)/2.d0
                  ezz=ezz+(vkl(1,3)**2+vkl(2,3)**2+vkl(3,3)**2)/2.d0
                  exy=exy+vkl(1,1)*vkl(1,2)+vkl(2,1)*vkl(2,2)+
     &                 vkl(3,1)*vkl(3,2)
                  exz=exz+vkl(1,1)*vkl(1,3)+vkl(2,1)*vkl(2,3)+
     &                 vkl(3,1)*vkl(3,3)
                  eyz=eyz+vkl(1,2)*vkl(1,3)+vkl(2,2)*vkl(2,3)+
     &                 vkl(3,2)*vkl(3,3)
!     
!     for frequency analysis or buckling with preload the
!     strains are calculated with respect to the deformed
!     configuration
!     
               elseif(iperturb(1).eq.1) then
                  exx=exx+vokl(1,1)*vkl(1,1)+vokl(2,1)*vkl(2,1)+
     &                 vokl(3,1)*vkl(3,1)
                  eyy=eyy+vokl(1,2)*vkl(1,2)+vokl(2,2)*vkl(2,2)+
     &                 vokl(3,2)*vkl(3,2)
                  ezz=ezz+vokl(1,3)*vkl(1,3)+vokl(2,3)*vkl(2,3)+
     &                 vokl(3,3)*vkl(3,3)
                  exy=exy+vokl(1,1)*vkl(1,2)+vokl(1,2)*vkl(1,1)+
     &                 vokl(2,1)*vkl(2,2)+vokl(2,2)*vkl(2,1)+
     &                 vokl(3,1)*vkl(3,2)+vokl(3,2)*vkl(3,1)
                  exz=exz+vokl(1,1)*vkl(1,3)+vokl(1,3)*vkl(1,1)+
     &                 vokl(2,1)*vkl(2,3)+vokl(2,3)*vkl(2,1)+
     &                 vokl(3,1)*vkl(3,3)+vokl(3,3)*vkl(3,1)
                  eyz=eyz+vokl(1,2)*vkl(1,3)+vokl(1,3)*vkl(1,2)+
     &                 vokl(2,2)*vkl(2,3)+vokl(2,3)*vkl(2,2)+
     &                 vokl(3,2)*vkl(3,3)+vokl(3,3)*vkl(3,2)
               endif
!     
!     storing the local strains
!     
               if(iperturb(1).ne.-1) then
                  eloc(1)=exx
                  eloc(2)=eyy
                  eloc(3)=ezz
                  eloc(4)=exy/2.d0
                  eloc(5)=exz/2.d0
                  eloc(6)=eyz/2.d0
               else
!     
!     linear iteration within a nonlinear increment:
!     
                  eloc(1)=vokl(1,1)+
     &                 (vokl(1,1)**2+vokl(2,1)**2+vokl(3,1)**2)/2.d0
                  eloc(2)=vokl(2,2)+
     &                 (vokl(1,2)**2+vokl(2,2)**2+vokl(3,2)**2)/2.d0
                  eloc(3)=vokl(3,3)+
     &                 (vokl(1,3)**2+vokl(2,3)**2+vokl(3,3)**2)/2.d0
                  eloc(4)=(vokl(1,2)+vokl(2,1)+vokl(1,1)*vokl(1,2)+
     &                 vokl(2,1)*vokl(2,2)+vokl(3,1)*vokl(3,2))/2.d0
                  eloc(5)=(vokl(1,3)+vokl(3,1)+vokl(1,1)*vokl(1,3)+
     &                 vokl(2,1)*vokl(2,3)+vokl(3,1)*vokl(3,3))/2.d0
                  eloc(6)=(vokl(2,3)+vokl(3,2)+vokl(1,2)*vokl(1,3)+
     &                 vokl(2,2)*vokl(2,3)+vokl(3,2)*vokl(3,3))/2.d0
               endif
!     
!     calculating the deformation gradient (needed to
!     convert the element stiffness matrix from spatial
!     coordinates to material coordinates
!     deformation plasticity)
!     
               if((kode.eq.-50).or.(kode.le.-100)) then
!     
!     calculating the deformation gradient
!     
c     Bernhardi start
                  xkl(1,1)=vkl(1,1)+1.0d0
                  xkl(2,2)=vkl(2,2)+1.0d0
                  xkl(3,3)=vkl(3,3)+1.0d0
c     Bernhardi end
                  xkl(1,2)=vkl(1,2)
                  xkl(1,3)=vkl(1,3)
                  xkl(2,3)=vkl(2,3)
                  xkl(2,1)=vkl(2,1)
                  xkl(3,1)=vkl(3,1)
                  xkl(3,2)=vkl(3,2)
!     
!     calculating the Jacobian
!     
                  vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))
     &                 -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1))
     &                 +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))
!     
!     inversion of the deformation gradient (only for
!     deformation plasticity)
!     
                  if(kode.eq.-50) then
!     
                     ckl(1,1)=(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))/vj
                     ckl(2,2)=(xkl(1,1)*xkl(3,3)-xkl(1,3)*xkl(3,1))/vj
                     ckl(3,3)=(xkl(1,1)*xkl(2,2)-xkl(1,2)*xkl(2,1))/vj
                     ckl(1,2)=(xkl(1,3)*xkl(3,2)-xkl(1,2)*xkl(3,3))/vj
                     ckl(1,3)=(xkl(1,2)*xkl(2,3)-xkl(2,2)*xkl(1,3))/vj
                     ckl(2,3)=(xkl(2,1)*xkl(1,3)-xkl(1,1)*xkl(2,3))/vj
                     ckl(2,1)=(xkl(3,1)*xkl(2,3)-xkl(2,1)*xkl(3,3))/vj
                     ckl(3,1)=(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))/vj
                     ckl(3,2)=(xkl(3,1)*xkl(1,2)-xkl(1,1)*xkl(3,2))/vj
!     
!     converting the Lagrangian strain into Eulerian
!     strain (only for deformation plasticity)
!     
                     cauchy=0
                     call str2mat(eloc,ckl,vj,cauchy)
                  endif
!     
               endif
!     
!     calculating fields for incremental plasticity
!     
               if(kode.le.-100) then
!     
!     calculating the deformation gradient at the
!     start of the increment
!     
!     calculating the displacement gradient at the
!     start of the increment
!     
                  do m2=1,3
                     do m3=1,3
                        vikl(m2,m3)=0.d0
                     enddo
                  enddo
!     
                  do m1=1,nope
                     do m2=1,3
                        do m3=1,3
                           vikl(m2,m3)=vikl(m2,m3)
     &                          +shp(m3,m1)*vini(m2,konl(m1))
                        enddo
                     enddo
                  enddo
!     
!     calculating the deformation gradient of the old
!     fields
!     
                  xikl(1,1)=vikl(1,1)+1
                  xikl(2,2)=vikl(2,2)+1.
                  xikl(3,3)=vikl(3,3)+1.
                  xikl(1,2)=vikl(1,2)
                  xikl(1,3)=vikl(1,3)
                  xikl(2,3)=vikl(2,3)
                  xikl(2,1)=vikl(2,1)
                  xikl(3,1)=vikl(3,1)
                  xikl(3,2)=vikl(3,2)
!     
!     calculating the Jacobian
!     
                  vij=xikl(1,1)*(xikl(2,2)*xikl(3,3)
     &                 -xikl(2,3)*xikl(3,2))
     &                 -xikl(1,2)*(xikl(2,1)*xikl(3,3)
     &                 -xikl(2,3)*xikl(3,1))
     &                 +xikl(1,3)*(xikl(2,1)*xikl(3,2)
     &                 -xikl(2,2)*xikl(3,1))
!     
!     stresses at the start of the increment
!     
                  do m1=1,6
                     stre(m1)=stiini(m1,jj,i)
                  enddo
!     
               endif
!     
!     prestress values
!     
               if(iprestr.eq.1) then
                  do kk=1,6
                     beta(kk)=-prestr(kk,jj,i)
                  enddo
               else
                  do kk=1,6
                     beta(kk)=0.d0
                  enddo
               endif
!     
               if(ithermal(1).ge.1) then
!     
!     calculating the temperature difference in
!     the integration point
!     
                  t0l=0.d0
                  t1l=0.d0
                  if(ithermal(1).eq.1) then
                     if((lakonl(4:5).eq.'8 ').or.
     &                    (lakonl(4:5).eq.'8I')) then
                        do i1=1,8
                           t0l=t0l+t0(konl(i1))/8.d0
                           t1l=t1l+t1(konl(i1))/8.d0
                        enddo
                     elseif((lakonl(4:6).eq.'20 ').or.
     &                       (lakonl(4:6).eq.'26 ')) then
                        nopered=20
                        call lintemp(t0,t1,konl,nopered,jj,t0l,t1l)
                     else
                        do i1=1,nope
                           t0l=t0l+shp(4,i1)*t0(konl(i1))
                           t1l=t1l+shp(4,i1)*t1(konl(i1))
                        enddo
                     endif
                  elseif(ithermal(1).ge.2) then
                     if((lakonl(4:5).eq.'8 ').or.
     &                    (lakonl(4:5).eq.'8I')) then
                        do i1=1,8
                           t0l=t0l+t0(konl(i1))/8.d0
                           t1l=t1l+vold(0,konl(i1))/8.d0
                        enddo
                     elseif((lakonl(4:6).eq.'20 ').or.
     &                       (lakonl(4:6).eq.'26 ')) then
                        nopered=20
                        call lintemp_th(t0,vold,konl,nopered,jj,t0l,t1l,
     &                         mi)
                     else
                        do i1=1,nope
                           t0l=t0l+shp(4,i1)*t0(konl(i1))
                           t1l=t1l+shp(4,i1)*vold(0,konl(i1))
                        enddo
                     endif
                  endif
                  tt=t1l-t0l
               endif
!     
!     calculating the coordinates of the integration point
!     for material orientation purposes (for cylindrical
!     coordinate systems)
!     
               if((iorien.gt.0).or.(kode.le.-100)) then
                  do j=1,3
                     pgauss(j)=0.d0
                     do i1=1,nope
                        pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1))
                     enddo
                  enddo
               endif
!     
!     material data; for linear elastic materials
!     this includes the calculation of the stiffness
!     matrix
!     
               istiff=0
!     
               call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,
     &              nalcon,imat,amat,iorien,pgauss,orab,ntmat_,
     &              elas,rho,i,ithermal,alzero,mattyp,t0l,t1l,ihyper,
     &              istiff,elconloc,eth,kode,plicon,nplicon,
     &              plkcon,nplkcon,npmat_,plconloc,mi(1),dtime,i,jj,
     &              xstiff,ncmat_)
!     
!     determining the mechanical strain
!     
               if(ithermal(1).ne.0) then
                  do m1=1,6
                     emec(m1)=eloc(m1)-eth(m1)
                  enddo
               else
                  do m1=1,6
                     emec(m1)=eloc(m1)
                  enddo
               endif
               if(kode.le.-100) then
                  do m1=1,6
                     emec0(m1)=emeini(m1,jj,i)
                  enddo
               endif
!     
!     subtracting the plastic initial strains
!     
               if(iprestr.eq.2) then
                  do m1=1,6
                     emec(m1)=emec(m1)-prestr(m1,jj,i)
                  enddo
               endif
!     
!     calculating the local stiffness and stress
!     
               call mechmodel(elconloc,elas,emec,kode,emec0,ithermal,
     &              icmd,beta,stre,xkl,ckl,vj,xikl,vij,
     &              plconloc,xstate,xstateini,ielas,
     &              amat,t1l,dtime,time,ttime,i,jj,nstate_,mi(1),
     &              iorien,pgauss,orab,eloc,mattyp,qa(3),istep,iinc,
     &              ipkon,nmethod,iperturb)
!     
               if(((nmethod.ne.4).or.(iperturb(1).ne.0)).and.
     &              (nmethod.ne.5)) then
                  do m1=1,21
                     xstiff(m1,jj,i)=elas(m1)
                  enddo
               endif
!     
               if(iperturb(1).eq.-1) then
!     
!     if the forced displacements were changed at
!     the start of a nonlinear step, the nodal
!     forces due do this displacements are 
!     calculated in a purely linear way, and
!     the first iteration is purely linear in order
!     to allow the displacements to redistribute
!     in a quasi-static way (only applies to
!     quasi-static analyses (*STATIC))
!     
                  eloc(1)=exx-vokl(1,1)
                  eloc(2)=eyy-vokl(2,2)
                  eloc(3)=ezz-vokl(3,3)
                  eloc(4)=exy-(vokl(1,2)+vokl(2,1))
                  eloc(5)=exz-(vokl(1,3)+vokl(3,1))
                  eloc(6)=eyz-(vokl(2,3)+vokl(3,2))
!     
                  if(mattyp.eq.1) then
                     e=elas(1)
                     un=elas(2)
                     um=e/(1.d0+un)
                     al=un*um/(1.d0-2.d0*un)
                     um=um/2.d0
                     am1=al*(eloc(1)+eloc(2)+eloc(3))
                     stre(1)=am1+2.d0*um*eloc(1)
                     stre(2)=am1+2.d0*um*eloc(2)
                     stre(3)=am1+2.d0*um*eloc(3)
                     stre(4)=um*eloc(4)
                     stre(5)=um*eloc(5)
                     stre(6)=um*eloc(6)
                  elseif(mattyp.eq.2) then
                     stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)
     &                    +eloc(3)*elas(4)
                     stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)
     &                    +eloc(3)*elas(5)
                     stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)
     &                    +eloc(3)*elas(6)
                     stre(4)=eloc(4)*elas(7)
                     stre(5)=eloc(5)*elas(8)
                     stre(6)=eloc(6)*elas(9)
                  elseif(mattyp.eq.3) then
                     stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)+
     &                    eloc(3)*elas(4)+eloc(4)*elas(7)+
     &                    eloc(5)*elas(11)+eloc(6)*elas(16)
                     stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)+
     &                    eloc(3)*elas(5)+eloc(4)*elas(8)+
     &                    eloc(5)*elas(12)+eloc(6)*elas(17)
                     stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)+
     &                    eloc(3)*elas(6)+eloc(4)*elas(9)+
     &                    eloc(5)*elas(13)+eloc(6)*elas(18)
                     stre(4)=eloc(1)*elas(7)+eloc(2)*elas(8)+
     &                    eloc(3)*elas(9)+eloc(4)*elas(10)+
     &                    eloc(5)*elas(14)+eloc(6)*elas(19)
                     stre(5)=eloc(1)*elas(11)+eloc(2)*elas(12)+
     &                    eloc(3)*elas(13)+eloc(4)*elas(14)+
     &                    eloc(5)*elas(15)+eloc(6)*elas(20)
                     stre(6)=eloc(1)*elas(16)+eloc(2)*elas(17)+
     &                    eloc(3)*elas(18)+eloc(4)*elas(19)+
     &                    eloc(5)*elas(20)+eloc(6)*elas(21)
                  endif
               endif
!     
!     calculation of the shape energy for the objective 
!     function: SHAPE ENERGY
!     
               ener(jj,i)=enerini(jj,i)+
     &              ((eloc(1)-eth(1)-emeini(1,jj,i))*
     &              (stre(1)+stiini(1,jj,i))+
     &              (eloc(2)-eth(2)-emeini(2,jj,i))*
     &              (stre(2)+stiini(2,jj,i))+
     &              (eloc(3)-eth(3)-emeini(3,jj,i))*
     &              (stre(3)+stiini(3,jj,i)))/2.d0+
     &         (eloc(4)-eth(4)-emeini(4,jj,i))*(stre(4)+stiini(4,jj,i))+
     &         (eloc(5)-eth(5)-emeini(5,jj,i))*(stre(5)+stiini(5,jj,i))+
     &         (eloc(6)-eth(6)-emeini(6,jj,i))*(stre(6)+stiini(6,jj,i))
               
               xenerel=xenerel+weight*xsj*ener(jj,i)
!     
!     end of loop over all integration points
            enddo
!     
!     summation of the objective function over all elements
!     
            if(idesvar.eq.0) then
               xener(i)=xenerel
               g0(iobject)=g0(iobject)+xenerel
            else
               dgdx(idesvar,iobject)=dgdx(idesvar,iobject)
     &              +(xenerel-xener(i))/distmin
            endif
!     
!     end of loop over all elements in thread
!     
         enddo
!     
!     
      return
      end
